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1. Introduction 

The Large Hadron Collider (LHC) is expected to start next year. The physics program of 
the two main experiments ATLAS and CMS aims at the discovery of the last missing particle 
predicted by the Standard Model - the Higgs boson - and at the search for signals of new physics 
beyond the Standard Model. Common to these searches is the fact that the signal events have to be 
digged out from a bulk of background events. The background events are due to Standard Model 
processes, mostly QCD processes which sometimes are accompanied by additional electro-weak 
bosons. Common to these searches is further the fact, that the final state is characterised by a high 
number of hadronic jets or identified particles. A rough theoretical understanding of scattering 
processes at hadron colliders is sketched in fig. |l|. The initial protons do not take part as a single 
entity in the scattering process, instead the basic constituents of the proton - partons like quarks 
and gluons - enter the hard scattering process. The probability of finding a specific parton inside 
the proton is described by parton distribution functions (pdf 's). The parton distribution functions 
are non-perturbative objects and therefore at present cannot be calculated from theory. However, 
they are universal and can be measured and extracted from one experiment and then used as input 
data for other experiments. In particular, the HERA experiment at DESY contributed significantly 
to our knowledge of the parton distribution functions. As a technical detail, the parton distribution 
functions depend on a scale. In simplified terms, they are measured at a scale Qq, but used as input 
data at a different scale 2i. In this context it is worth to note that the variation with the scale of the 
parton distribution functions can be calculated within perturbation theory. 

The hard scattering process can be formulated entirely in terms of the fundamental fields 
(quarks, gluons, ...) of the Standard Model. It is calculable in perturbation theory. Attached to 
the hard scattering process are parton showers, where the partons of the hard scattering event radi- 
ate off additional coUinear or soft partons. From this additional radiation originates the observed 
hadronic jets, i.e. bunches of particles moving in the same direction. The parton shower increases 
significantly the number of partons in the event. In principle, the parton shower is governed to a 
large extend by perturbation theory. But due to our limited computational abilities we are forced 
to replace the full matrix elements by approximations. These approximations are based on the 
observation, that matrix elements are enhanced in soft and collinear regions. 

After the parton shower, the partons convert to hadrons. Like for the parton distribution func- 
tions this is non-perturbative physics. In practise, models like the string model or the cluster model 
are used to mimic the formation of hadrons. If unstable hadrons are formed, they subsequently 
decay. 

This simplified picture of an event has to be completed with additional complications due to 
the underlying event, multiple interactions and pile-ups. As the hard scattering breaks up the initial 
protons, the proton remnants are not colour-neutral. The underlying event describes the interactions 
of the proton remnants. Usually the underlying event will produce activity in the detector along 
the beam pipe. However, it can happen that more than one pair of partons of the initial protons 
undergo a hard scattering. This goes under the name of multiple interactions. Finally, there is 
the possibility that more than one proton-proton scattering takes place in a single bunch crossing. 
These are pile-up events. 

As can be seen, the theoretical description of a single event is rather involved and requires 
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Figure 1: A schematic and simplified description of an event in hadron-hadron collisions. 

in particular knowledge on non-perturbative physics. As this knowledge is not available, we have 
to content ourselves to models and approximations. This will limit the accuracy and precision of 
theoretical predictions. Fortunately, for a specific class of observables we can do better: This is the 
class of infrared-safe observables. Infrared safeness implies that the value of an observable does 
not change, if additional soft or coUinear particles are added to an event: 

lim 0„+r{pi,...,p„+r) = 0„{p[,...,p'„) (1.1) 

r partons unresolved 

Here, 0„ (pi, ...,p„) denotes the value of the observable for an event with particles with four- 
momenta pi, p„. Infrared-safe observables depend only mildly on showering and hadronisation 
and can be calculated reliably in perturbation theory. At hadron colliders we need of course the 
additional non-perturbative information on the parton distribution functions, but these quantities 
have been measured and are available. The master formula for the calculation of an observable is 
given by 

[d<l>n-20{pu.:,Pn) Kl'- (1-2) 

The various ingredients of this formula are: The parton distribution functions, i.e. the probability of 
finding a parton i with momentum fraction x inside the parent hadron h are denoted by fi{x). 2K{s) 
is the flux factor, equal to two times the center-of-mass energy squared of the two incoming partons. 
l/(2/i + l)/(2/2 + l)/ni/n2 corresponds to an averaging over the spins and colour degrees of 
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Figure 2: Modelling of jets in perturbation theory. At leading-order a jet is modelled by a single parton, at 
next-to-leading order either by one or two partons. At next-to-next-to-leading order a jet is modelled by one, 
two or three partons. 

freedom for the initial particles. The second sum is over the number of final-state particles, the 
integral is over the phase space of {n — 2) final-state particles. The most important ingredient of 
this formula is the matrix element squared for « particles, (« — 2) in the final state, two in 
the initial state. As can be seen from this formula, each event is weighted by the matrix element 
squared. 

At high energies the strong coupling is small and the matrix element can be calculated in 
perturbation theory. For an observable, whose leading-order contribution in perturbation theory is 
given by an n-parton amplitude, the following expansions are relevant for the calculation of the 
next-to-leading order (NLO) and next-to-next-to-leading order (NNLO) predictions: 
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(1.3) 



In this formulae 



A') 



denotes an amplitude with n external particles and / loops. At leading- 
order only the Bom amplitude j^n^^ contributes. At next-to-leading order we have in addition the 
contributions from the one-loop ampUtude ^n^^ and the single emission contribution, which is 
given by the Bom ampUtudes <<2^^+i with one additional parton. These two contribution are of the 
same order with respect to the power counting of the coupling. As far as the phase space for the 
final-state particles is concerned, these two contributions live on different phase spaces of different 
dimensions, since the single emission contribution has one additional particle in the final state. At 
next-to-next-to-leading order we have to include in addition the contributions from the two-loop 
ampUtude and the one-loop ampUtude squared, the one-loop ampUtudes with a single additional 
emission and the Bom amplitudes with two additional emissions. Up to next-to-next-to-leading 
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Figure 3: Scale-dependence of the cross-section for the process pp ff +jet. The leading-order pre- 
diction shows a strong scale-dependence. The next-to-leading order prediction reduces significantly the 
scale-dependence. The plot is taken from [1]. 

order we include the radiation of up to two additional particles. Therefore jets are modelled by 
one, two or three partons. This is shown in fig. ^ 

Why are the computations of higher-order corrections necessary ? The answer is that we want 
to achieve a certain precision. This is illustrated by a simple example. Let us consider a pure QCD 
process with three hard partons in the final state. The leading-order prediction is proportional to 
a^. It is a well-known fact that the numerical value of the strong coupling depends on an arbitrary 
scale jXren, at leading order the formula reads 
471 2 

a, = —, j3o = 1 1 - -Nf, A « 165 MeV for Nf = 5 hght flavours. (1.4) 

jSoln^ 3 

As the choice of the scale is arbitrary, and since enters the leading-order prediction to the third 
power, this introduces a strong scale-dependence of the theoretical prediction and therefore a large 
uncertainty. The inclusion of higher-order corrections reduces this uncertainty. This is shown for 
the example of tt + jet production at the LHC in fig. ^. This is an example where the leading-order 
prediction is proportional to a^. Shown in fig. ^ is the variation of the theoretical prediction with 
the scale jJ.. Apart from the scale entering the formula for the strong coupling, which is called the 
renormalisation scale there is in addition also a scale, at which the parton distribution functions are 
evaluated. The latter is called the factorisation scale. In the plot both scale have been identified, 

IJ- = IJ-ren = IJ-fact ■ (1-5) 

As can be seen, the inclusion of the next-to-leading order correction reduces significantly the scale- 
dependence and clearly motivates the quest for the computation of higher-order corrections. 

What are the objectives for LHC physics ? Let us now look more closely where precision 
calculations are needed for LHC physics. There are three important points to mention: 

• In order to predict absolute rates to a good precision, we have to know the input parameters 
of a theoretical calculation to a high precision. The most important input parameters are the 
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process relevant for 



1. 


PP - 


VV + jet 


ttH, new physics 


2. 


PP - 


+ 2 jets 


Higgs production by vector boson fusion 


3. 


PP- 


-y ttbb 


ttH 


4. 


PP- 


-y « + 2 jets 


ttH 


5. 


PP - 


-^VVbb 


VBF ^H^VV, ttH, new physics 


6. 


PP - 


VV + 2jets 


VBF -^H^VV 


7. 


PP - 


V + 3jets 


various new physics signatures 


8. 


PP - 


-^VVV 


SUSY 



Table 1: The experimenter's wish Hst at Les Houches 2005 [2]. V denotes either one of the electro-weak 
bosons: photon, Z-boson or W-boson. 

strong coupling and the parton distribution functions. These parameters can be extracted 
from other experiments, for example the strong coupling can be extracted from the process 
e+e^ —>■ 3 jets at LEP, the pdf 's are measured at HERA. The experimental precision has 
to be matched with the corresponding precision of the theoretical calculation, requiring the 
inclusion of the NNLO corrections. This in turn requires a calculation for e^e^ —>■ 3 jets 
at NNLO and the computation of the three-loop splitting functions for the evolution of the 
pdf's. 

• Standard candles at the LHC: For a few standard hard pp processes like the production 
of W, Z, top, Higgs or jets we would like to have a precise prediction from theory, again 
requiring a NNLO computation. Higgs production is of clear interest for the search of the 
last undiscovered particle predicted by the Standard Model. The other standard processes are 
useful to measure fundamental quantities like the W- or top-mass to a better precision or can 
be used by the experimentalists to understand their detector. 

• Finally one would like to have reliable predictions for multi-particle final states that occur at 
high rates and form background to new physics. This implies next-to-leading order calcula- 
tions. Examples for such processes have been collected in the experimenter's wish list at Les 
Houches workshop in 2005 and are shown in table [l[ 

What is the state-of-the-art for multi-leg processes ? It is clear that the complexity of a calcu- 
lation increases with the order in perturbation theory as well as with the number of external legs. 
Therefore we expect that if we go up in the order of perturbation theory the available calculations 
are restricted to fewer final state particles. Let us first consider leading order calculations. At this 
order there are techniques for the automated calculation of the (leading order) matrix elements 
[3-10] as well as techniques for the efficient integration over phase space [11, 12]. Several com- 
puter programs like Madgraph/Madevent [13-15], Sherpa/Amegic-i~i- [16], Helac/Phegas [17, 18], 
Comphep [19], Grace [20] or Alpgen [21] are available, which implement these techniques and can 
be used to obtain leading order predictions for processes with a rather high number of final state 
particles. 
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n 


2 


3 


4 


5 


6 


7 


8 


diagrams 


4 


25 


220 


2485 


34300 


559405 


10525900 



Table 2: Number of Feynman diagrams contributing to gg ng at tree level. 



If we now consider the next order in perturbation theory, we observe that there many NLO 
calculation for 2 ^ 2 processes at hadron colliders, but only a few for 2 ^ 3 processes. Fully 
differential numerical programs for 2 — > 3 processes exist for example for pp — > 3 jets [22-24], 
pp V + 2 jets [25-29], pp ttH [30, 31], pp ^ // + 2jets [32-34], pp «+jet [1] and 
pp ^ ZZZ [35]. Of comparable complexity is the NLO calculation for the production of two vector 
bosons in vector-boson-fusion [36-38]. On the other side there aren't at present any NLO programs 
for the LHC with more than 3 hard particles in the final state. NLO programs with four particles in 
the final state are available for electron-positron annihilation for e+e^ 4 fermions [39, 40] and 
e^e^ —>■ 4 jets [41-44]. As can be seen, for NLO computations our capabilities are already much 
more restricted with respect to the number of final state particles. For processes with only two 
or three final state particles tools which help with the automatisation of the computation like the 
combination FeynArts, FormCalc, Looptools [45-48] or the Grace package [49] are available. 

Finally, we consider NNLO calculations. Here, fully differential predictions for the LHC are 
available only for a few selected 2^1 and 2^2 processes, like Drell-Yan [50], -production [51] 
or Higgs production [52-58]. For electron-positron annihilation NNLO predictions are available 
for e^e^ —>■ 2 jets [59-61] and the thrust distribution [62]. 

What are the bottle-necks ? As can be seen from this summary of the state-of-the-art, higher- 
order corrections are limited to processes with not too many particles in the final state. Let us now 
look into the difficulties, which prohibit a straight-forward automatisation of the computation of 
higher order corrections for processes with many particles in the final state. These difficulties can 
be grouped into three categories: 

• Length: Perturbative calculations lead to expressions with a huge number of terms. 

• Integrals: At one-loop and beyond, the occurring integrals cannot be simply looked up in an 
integral table. 

• Divergences: At NLO and beyond, infrared divergences occur in intermediate stages, if 
massless particles are involved. 

The first complication - lengthy expressions - affects already leading-order computations. I will 
discuss methods to handle this problem in section ^. The second and the third item occur for the 
first time in NLO computations. They are discussed in section ^ and section Q respectively. 

2. Managing lengthy expressions 

It is a well-known fact that the complexity of a calculation based on Feynman diagrams growth 
factorially with the number of external particles. As an example consider the all-gluon amplitude 
at tree level. Table |^ shows the number of diagrams contributing to the process gg ng at tree 
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level. All diagrams involve three- and four-gluon vertices. The Feynman rules for these vertices 
blow-up the corresponding expressions even further: 




= gf''' [ih - h)^gvx + {h - h)vgxii + {h - h)xg^v] , 




= -ig^ [f"''f"' [g^^xgyp -g^pgvx) +r'f^' {g^.g^p -g^pgxv) 

+f'"r'{s^^gxp-g,xgvp)] 

(2.1) 

For the computation of observables we have to square the amplitude and to sum over all spins or 
helicities. For gluons we replace the polarisation sum by 

£^[k,A)ev[k,X) = -g^v + ^ ^ Jknf' ^ 

In this formula, is an arbitrary four- vector. If the amplitude consists of 0{N) terms, squaring 
the amplitude will produce 0{N^) terms. From this example it is clear, that the number of terms 
in intermediate expressions of a calculation based on Feynman diagrams growth dramatically with 
the number of external legs. As already mentioned, this problem occurs already at tree level. 

2.1 Computer algebra 

As calculations for multi-particle final states tend to involve lengthy intermediate expressions, 
computer algebra has become an essential tool. In fact, particle physics is and has been a driving 
force for the development of computer algebra systems. Quite a few computer algebra systems 
have their roots within the high energy physics community or strong links with them: REDUCE, 
SCHOONSHIP, MATHEMATICA, FORM [63] or GiNaC [64, 65], to name only a few. In most 
cases the requirements on a computer-algebra system for computer-intensive symbolic calculations 
in particle physics can be summarised by: 

• The computer algebra system has to provide basic operations like addition, multiplication, 
sorting, etc.. 

• Specialised code for the solution of a particular problem is usually written by the user. The 
computer algebra system has to provide a convenient programming language. 

• There is actually no need for a system which knows "more" than the user. 

The most widely used computer algebra systems in the community are MATHEMATICA, MAPLE, 
REDUCE, FORM and GiNaC. The first three are commercial programs, FORM and GiNAC are 
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non-commercial and freely available . Below there are two small example programs in FORM and 
GiNaC for the calculation of 

Tr ^i^2^3^^4 = 4 (pi • P2) {p3 ■p4)-4{pi-p3) {P2 • P4) + 4 (pi • pa) {P2 -ps)- (2.3) 
The example in FORM reads: 

* Example program for FORM 
V pl,p2,p3,p4; 

L res = g_(l,pl) *g_(l,p2) *g_(l,p3) *g_(l,p4) ; 
trace4 , 1 ; 
print ; 
. end 

GiNaC is a C++ Ubrary, which provides capabiUties for symbolic calculations within the C++ 
programming language. The corresponding example in GiNaC reads: 

#include <iostream> 
#include <ginac/ginac. h> 

using namespace std; 
using namespace GiNaC; 

int main ( ) 
{ 

varidx mu ( symbol ( "mu" ), 4 ) , nu ( symbol ( "nu" ) ^ 4 ) , 

rho ( symbol ( " rho" ) , 4 ) , sigma ( symbol ( " sigma" ) , 4 ) ; 

symbol pi ("pi"), p2("p2"), p3("p3"), p4("p4"); 

ex res = dirac_gamma (mu, 1 ) *dirac_gamma (nu, 1 ) 
*dirac_gamma (rho, 1) *dirac_gamma (sigma, 1) 

*indexed (pi , mu . toggle_variance ( ) ) *indexed (p2 , nu . toggle_variance ( ) ) 

* indexed (p3, rho . toggle_var lance ( ) ) * indexed (p4 , sigma . toggle_variance ( ) ) ; 

res = dirac_trace ( res , 1 ) ; 

res = res . expand ; 

res = res . simplif y_indexed ( ) ; 

cout << res << endl; 

} 

Computer algebra is an essential tool, but a brute force application alone will still produce 
lengthy expressions, which are slow and potentially unstable, when evaluated numerically. 

^FORM is available at http://www.nikhef.nirfomi, GiNaC is available at http://www.ginac.de. 
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2.2 Quantum number management 

In order to keep the size of intermediate expressions under control, a divide-and-conquer strat- 
egy has been proven useful: One divides the quantity to be calculated into smaller pieces and 
calculates the small pieces separately. This approach is also called "quantum number manage- 
ment". One first observes that it is not necessary to square the amplitude and sum over the spins 
and helicities analytically. It is sufficient to do this numerically. This avoids obtaining 0{N^) terms 
from an expression with 0{N) terms. The individual amplitudes have to be calculated in a helicity 



or spin basis. This is discussed in section [2.2. 1| . The second observation is related to the fact, that 
individual helicity amplitudes can be decomposed into smaller gauge-invariant pieces, called par- 
tial amplitudes. This is discussed in section 2.2.2| . These partial amplitudes can be calculated for a 



given helicity configuration without reference to Feynman diagrams, as discussed in section 2.2.3 . 
Finally, for specific helicity configurations, compact formula are known, which are summarised 
in section 



2.2.4. Supersymmetric identities provide relations between amplitudes with different 
particle contents, therefore only some of them need to be calculated. This is discussed in section 
12j. 



2.2.1 Helicity amplitudes 

The computation of helicity amplitudes requires the choice of a helicity (or spin) basis. This 
is straightforward for massless fermions. The two-component Weyl spinors provide a convenient 
basis: 

\p±) = ^{l±r5)u{p). (2.4) 

In the literature there are different notations for Weyl spinors. Apart from the bra-ket-notation 
there is the notation with dotted and un-dotted indices: The relation between the two notations is 
the following: 

\p+)=Pb, {p + \=Pa, 

\P-)=P^ {P-\=P^- (2.5) 

Spinor products are defined as 

{pi) = {p - \q+) , [pq\ = {p + \q-), (2-6) 

and take value in the complex numbers. It was a major break-through, when it was realised that 
also gluon polarisation vectors can be expressed in terms of two-component Weyl spinors [66-71]. 
The polarisation vectors of external gluons can be chosen as 

where k is the momentum of the gluon and q is an arbitrary light-like reference momentum. The 
dependence on the arbitrary reference momentum q will drop out in gauge invariant quantities. The 
polarisation sum is that of an light-like axial gauge: 

£ e',{k,q) [e'Ak.q))* = -^,v + ^^^^^^^i^. (2.8) 

A— i ^ 
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Changing the reference momentum will give a term proportional to the momentum of the gluon: 

e;(t,„)-e;(*,..) = V2J|l|l^.,. (2.9) 



For massive fermions we can take the spinors as [72, 73] 



Here is a light-Uke four vector obtained through 

p^=p-^q. (2.11) 
2pq 

q denotes again an arbitrary hght-like reference momentum and is related to the quantisation axis 
of the spin for the massive fermion. In contrast to the gluon case individual ampUtudes with label 
+ or — will depend on the choice of the reference momentum q. 

2.2.2 Colour decomposition 

Throughout this article I use the normalisation 

TrTT* = ig"^ (2.12) 

for the colour matrices. Amplitudes in QCD may be decomposed into group-theoretical factors 
(carrying the colour structures) multipUed by kinematic functions called partial amplitudes [74—78]. 
These partial ampUtudes do not contain any colour information and are gauge-invariant objects. In 
the pure gluonic case tree level ampUtudes with n external gluons may be written in the form 

K{h2,...,n) =g"-^ £ 2Tr(r«-W...r«''('.))A„((7(l),...,(7(n)), (2.13) 

OeSn/Zr, 

where the sum is over all non-cyclic permutations of the external gluon legs. The quantities 
A„(a(l), a(?i)), called the partial amplitudes, contain the kinematic information. They are 
colour-ordered, e.g. only diagrams with a particular cycUc ordering of the gluons contribute. The 
colour decomposition is obtained by replacing the structure constants f"'"^ by 



if be ^ 2 Tr (t"T^T'^ -Tr (t''T''T'^ 



(2.14) 



which follows from [r°,r*] = if^^T^. The resulting traces and strings of colour matrices can be 
further simpUfied with the help of the Fierz identity : 

TuTSi = \{^^n^jk-]^5iA^- (2.15) 
The colour decomposition for a tree ampUtude with a pair of quarks is 

K+2{q,lX...,n,q) = g« £ (r'^W ...r-W),^^. A„+2(^, a(l), a(2), a(n),^-). (2.16) 

Sn 
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where the sum is over all permutations of the gluon legs. In general, the colour factors are combi- 
nations of open 

(T"'...!""), ■ (2.17) 

and closed strings 

Trl^r^'.-.r^") (2.18) 

of colour matrices. These building blocks form a basis in colour space. The choice of the basis for 
the colour structures is not unique, and several proposals for bases can be found in the literature 
[79-81]. A second useful basis is the colour-flow basis: This basis is obtained by replacing every 
contraction over an index in the adjoint representation by two contractions over indices i and j in 
the fundamental representation: 

ya^a ^ yagab^b ^ ya (^T^^.T^^j = (VlTf'jV"^ (VlTfiE^^ . (2.19) 

In this representation the colour decomposition of the pure gluon ampUtude is given by 

<(l,2,...,n) = £ 5,,^,,^5^,,^...5,,„,,A„((Ji,...,(7„). (2.20) 

vv2y ^^s„iz, 

2.2.3 Off-shell recurrence relations 

Berends-Giele type recurrence relations [3, 82] build partial ampUtudes from smaller building 
blocks, usually called colour-ordered off-shell currents. Off-shell currents are objects with n on- 

shell legs and one additional leg off-shell. Momentum conservation is satisfied. It should be noted 
that off-shell currents are not gauge-invariant objects. Recurrence relations relate off-shell currents 
with n legs to off-shell currents with fewer legs. The recursion starts with n = l: 

J^'ih) = e^'ikuq). (2.21) 

is the polarisation vector of the gluon and q an arbitrary light-like reference momentum. The 
recursive relation states that in the pure-gluon off-shell current a gluon couples to other gluons only 
via the three- or four-gluon vertices : 



J^{k^\ ...,k^) — 2 



"f y3'^^^(-*ri,„,*:ij,*:,-+i,„)yv(fcf',...,40-^p(^>r/'---'^^) 



n—2 n—\ 



+ 22 H ^ ■,----,k/)Jp{kj+i-i----iki^)Ja{ki!^i,---,k^) 

j=U=j+l 



,(2.22) 



where 

Kij = ki + ki+i + ... + kj (2.23) 
and ^3 and V4 are the colour-ordered three-gluon and four-gluon vertices 
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Figure 4: Off-shell recurrence relation: In an off-shell current particle n + 1 is kept off-shell. This allows to 
express an off-shell current with n on-shell legs in terms of currents with fewer legs. 

Vr'"" = /(Ig^Pg^^ -g'' -g'^ V) • (2.24) 
The recurrence relation is shown pictorially in fig. ^ The gluon current is conserved: 

(tk^^Jn=^- (2-25) 

From an off-shell current one easily recovers the on-shell amplitude by removing the extra propa- 
gator, taking the leg (« + 1) on-shell and contracting with the appropriate polarisation vector. 

2.2.4 Parke-Taylor formulae 

The partial amplitudes have for specific helicity combinations remarkably simple analytic for- 
mula or vanish altogether. For the all-gluon tree amplitude one finds 

A„(l+,2+,...,n+) = 0, 
A„(l+,2+,...,r,-,«+) = 0, 

/ ^\ n-2 I jJc)^ 

A„{\+X,-,r,-,k-,...,n+) = i[V2) J^-^- (2-26) 

The partial amplitudes where all gluons have positive helicities, or where all gluons except one 
have positive helicities vanish. The first non-vanishing result is obtained for the n-gluon amplitude 
with n — 2 gluons of positive helicity and 2 gluons of negative helicity. It is given by a remarkable 
simple formula. Note that this formula holds for all n. An amplitude with n — 2 gluons of positive 
helicity and 2 gluons of negative helicity is called a maximal-helicity violating amplitude (MHV 
amplitude). Obviously, we find similar formula if we exchange all positive and negative helicities: 

A„(l-,2-, ...,?!-) = 0, 
A„(r,2-,...,7+, ...,?!-) = 0, 

A,^ira-,...j\...,k-,...,n-)=i{V2y-'^-^^ (2.27) 

These formulae have been conjectured by Parke and Taylor [83] and have been proven by Berends 
and Giele [3]. 
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2.2.5 Supersymmetric relations 

After removing the colour factors, QCD at tiee-level may be viewed as an effective supersym- 
metric theory [84-88], where the quarks and the gluons form a super-multiplet (a A'^ = 1 vector 
super-multiple). Let us denote by 

QsusY = e(q^QA+qAQ^) (2.28) 

the SUSY generators contracted with two-component Weyl spinors and qj^ and multiplied by 
a Grassmann number 0. In an unbroken supersymmetric theory, the supercharge annihilates the 
vacuum, and therefore 

n 

(0|[esusY,<J>i<J>2...<J>J|0) = X(0|*fi-[esusY,<J>,-]...<J>«|0) =0 (2.29) 

!=1 

where the field denotes either a gauge boson § or a fermion A. The commutators are given by 

[esusY,g^(^)] =r±(^,^)A±(^), [Q,^sY,A'^{k)]=r^{k,q)g^{k), (2.30) 

with 

r^{k,q) = d{q±\kT). (2.31) 

Let us now consider 



= ^^{Pl,q)An{gt,gt,--;gJ,--;gt-l'Sn)-'^^iPj,(l)MK'g2'---'-^J^--^8n-l'gn) 

-r-{pn,q)An{A+,g+,...,gJ,...,g+_„A-). (2.32) 

If one further sets the reference momentum equal to ^ = pj and uses the expression for the max- 
imally helicity violating gluon amplitudes one obtains the expression for an amplitude with a pair 
of quarks: 

An{qr,gi,.:,gJ,:;gUqn) = ' (^) (I^^^^ " ^^.33) 



2.3 New developments: Twistor methods 



In the previous section we saw that for the all-gluon tree amplitude there are remarkable simple 
formula if almost all gluons have one helicity and not more than two gluons have the opposite 
helicity. Of course we are also interested in the case, where more than two gluons have the opposite 



helicity. MHV vertices, discussed in section 2.3.1, provide an answer to this question and tell 



us how the complexity increases with the number of opposite helicity gluons. In addition this 
construction led to new recursion relation, which no longer require that one external leg is kept 
off-shell. The building blocks of these on-shell recursion relations are gauge-invariant amplitudes. 



This is discussed in section 2.3.2 
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Figure 5: MHV diagrams contributing to the tree-level six-gluon amplitude A6(l ,2 ,3 ,4+,5+,6+). 



2.3.1 MHV vertices 

As an alternative to usual Feynman graphs, tree amplitudes in Yang-Mills theory can be con- 
structed from tree graphs in which the vertices are tree level MHV scattering amplitudes, continued 
off shell in a particular fashion [89]. The basic building blocks are the MHV amplitudes, which 
serve as new vertices: 

Vnil+,...,j-,...,k-,...,n+) = i{'^y' \y2)''^'^l) - ^^-^^^ 

Each MHV vertex has exactly two lines carrying negative helicity and at least one line carrying 
positive helicity. Each internal line has a positive helicity label on one side and a negative helicity 
label on the other side. The propagator for each internal line is the propagator of a scalar particle: 

^ (2.35) 



The expression (2.34) for the MHV vertices involves spinors corresponding to massless on-shell 



momenta kj = 0. Therefore we have state what this light-like four-vector should be for every 
internal line meeting a MHV vertex. As in eq. (2. 1 1 ) the light-like four- vector can be taken as [90] 

k^ = k-:^q, (2.36) 
2k ■ q 

where k is the momentum flowing through the internal Une and g is a fixed light-like reference 
momentum. Let us now consider an example. The amplitude A6(l^,2^,3^,4+,5+,6+) has three 
gluons of positive helicity and three gluons of negative helicity and is one of the first non-trivial 
amplitudes, which are non-zero and which are not MHV amplitudes. Fig. ^ shows the six MHV 
diagrams contributing to this amplitude. The first diagram yields 



4+ 5+ 6+ 




<: 



(i2)(2(-^y)((-4)i) 



12 



V2 



{3k\^) 



{34){45){56){6k\,){k\,3) 



(2.37) 
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Figure 6: Diagrams contributing to the tree-level six-gluon amplitude A6(l ,2 ,3 ,4+,5+,6+) in the on- 
shell recursive approach. The vertices are on-shell amplitudes. 

and similar expressions are obtained for the five other diagrams. In this expression, ku = ki +k2 
is the momentum flowing through the internal line and k\2 is the projection onto a light-like four- 
vector as in eq. (2.36). We recall from table |2| that a brute force approach would require the 



calculation of 220 Feynman diagrams. Restricting ourselves to a partial amplitude with a fixed 
cychc order reduces this number to 36 diagrams. In the approach based on MHV vertices there are 
only six diagrams. In the next subsection we will discuss a method which reduces the number of 
diagrams even further. 

2.3.2 On-shell recursion relations 

Britto, Cachazo and Feng [91] gave a recursion relation for the calculation of the «-gluon 
amplitude which involves only on-shell amplitudes. To describe this method it is best not to view 
the partial amplitude A„ as a function of the four-momenta , but to replace each four- vector by a 
pair of two-component Weyl spinors. In detail this is done as follows: Each four-vector has a 
bi-spinor representation, given by 

^AB = ^tl'^AB-' ^ti = 2^AB^^^- (2.38) 

For light-like vectors this bi-spinor representation factorises into a dyad of Weyl spinors: 

k^k^ =0 ^ k^B = kAkg. (2.39) 



The equations (2.38) and ( 2.39| ) allow us to convert any light-like four-vector into a dyad of Weyl 



spinors and vice versa. Therefore the partial amplitude A„, being originally a function of the mo- 
menta kj and helicities Ay, can equally be viewed as a function of the Weyl spinors k^^, k^^ and the 
helicities Ay: 

) = ^nl^A'^lj-^l' ■••'^A'^^'-^n)- (2.40) 

Note that for an arbitrary pair of Weyl spinors, the corresponding four-vector will in general be 
complex- valued. If (A„, Ai) / (+, — ) we have the following recurrence relation: 

An{klk\,Xu-rAXBX) = (2.41) 

^ ^ Ay Ai,^^,^^, A2, ,Xj-i,iKA,iK^,—X 

7=3 A=± 

I . / . . ; . ; . ,„_i ,„_i 



■ „2 ^«-7'+2 ( -^Ai-^B) A,/c^,A;^, Ay, ,kg ,Xn-l,kA,kg,Xn 
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method 


diagrams 


brute force approach 


220 


colour-ordered amplitudes 


36 


MHV vertices 


6 


on-shell recursion 


2 



Table 3: The number of diagrams contributing to the colour-ordered six-gluon amplitude 

A6(l^,2^,3^,4+,5+,6+) using various methods. 



If Ai) = (+, — ) we can always cyclic permute the arguments, such that Ai) 7^ (+, — ). This 
is possible, since on-shell amplitudes, where all gluons have the same helicity, vanish. In eq. ( 2.41 ) 
the shifted spinors k^,k"j^, Ka and are given by 

K 

kA = k^ — zk^ , Ka 



k^ = k^ + zkl, K,= 'l^f (2.42) 



where 



^Afi=IW' ^iVi=det^A5, and z = ^--i^±-^. (2.43) 

Let us again consider as an example the amplitude A6(l^,2^,3^,4+,5+,6+). The diagrams con- 
tributing in the on-shell approach are shown in fig. @. One obtains as a result for the amplitude 



A6(r,2-,3-,4+,5+,6+) = 

(6 + |l+2|3+)3 , (4 + |5 + 6|l+)- 



4/ 



■ + ■ 



(2.44) 



[61] [12] (34) (45)^126 (2 + 1 1 + 6|5+) [23] [34] (56) (61)si56 (2 + 1 1 + 6|5+) _ 

Note that there only two diagrams, which need to be calculated. Table ^ shows a comparison 
of the number of diagrams contributing to the colour-ordered six-gluon amplitude in the various 
approaches. The performance of a numerical implementation of these new methods have been 
investigated in [92,93]. 

2.3.3 Proof of the on-shell recursion relations 

For the proof [94-99] of the on-shell recursion relation one considers the function 

A{z) =An{k\,klXi,...,ki;kl,X„) (2.45) 

of one variable z, where the z-dependence enters through 

k^^=k\-zkl, kl = kl + zkl. (2.46) 

Note that for all values of z, this is an on-shell amplitude. However, the four-momenta of par- 
ticles 1 and n are in general complex. The function A{z) is a rational function of z, which has 
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only simple poles in z- This follows from the Feynman rules and the factorisation properties of 
amplitudes. Therefore, if A(z) vanishes for z —>-°°, A{z) is given by Cauchy's theorem as the sum 
over its residues. This is just the right hand side of the recursion relation. The essential ingre- 
dient for the proof is the vanishing of A(z) at z ^ °°. If = (+,— ) it can be shown that 
each individual Feynman diagram vanishes for z ^ °°- Consider the flow of the z-dependence in 
a particular diagram The most dangerous contribution comes from a path, where all vertices are 
three-gluon-vertices. For a path made of n propagators we have n + l vertices and the product of 
propagators and vertices behaves therefore like z for large z- This statement remains true for a path 
containing only one vertex and no propagators. The polarisation vectors for the helicity combina- 
tion (Ai,A„) = (+,— ) contribute a factor 1/z^, therefore the complete diagram behaves like 1/z 
and vanishes therefore for z^ °°. 



3. Calculating loop amplitudes 

The second bottle-neck for higher-order computations are loop integrals. I first discuss one- 



loop integrals in section 3.1. These are relevant for multi-leg NLO calculations, like the processes 



listed in table In section 3.2 I discuss techniques for two-loop amplitudes and beyond. 



3.1 Automated computation of one-loop amplitudes 

The simplest, but most important loop integrals are the one-loop integrals. We have a good 
understanding of these integrals and I will present the main results in this section. An important 
result is that any scalar integral with more than four external legs can be reduced to scalar inte- 
grals with no more than 4 external legs. Therefore the set of basic one-loop integrals is rather 
limited. I will discuss this reduction in section p.l.l| . For one-loop tensor integrals we can use the 
Passarino-Veltman method, which reduces any tensor integral to a combination of scalar integrals. 
Improvements of the Passarino-Veltman algorithm are discussed in section [3.1.2[ Section [3.1.3| is 
devoted to methods, which avoid Feynman diagrams. 

3.1.1 Reduction to integrals with no more than four external legs 

In this section we discuss the reduction of scalar integrals with more than four external legs to 
a basic set of scalar one-, two-, three- and four-point functions. It is a long known fact, that higher 
point scalar integrals can be expressed in terms of this basic set [100, 101], however the practical 
implementation within dimensional regularisation was only worked out recently [102-106]. The 
one-loop n-point functions with « > 5 are always UV-finite, but they may have IR-divergences. Let 
us first assume that there are no IR-divergences. Then the integral is finite and can be performed in 
four dimensions. In a space of four dimensions we can have no more than four linearly independent 
vectors, therefore it comes to no surprise that in an one-loop integral with five or more propagators, 
one propagator can be expressed through the remaining ones. This is the basic idea for the reduction 
of the higher point scalar integrals. With slight modifications it can be generalised to dimensional 
regularisation. I will discuss the method for massless one-loop integrals 

/„ = e-y^^i^-{-ir / —-^ ^ ^. (3.1) 

J iK2 k^{k- piY..\k- Pi - ...pn-\Y 
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With the notation 



one can associate two matrices S and G to the integral in eq. (3.1). The entries of the nxn kine- 
matical matrix S are given by 

Sij = {qi-qjf, (3.3) 

and the entries of the (« — 1) x (« — 1) Gram matrix are defined by 

Gij = Iqiqj. (3.4) 

For the reduction one distinguishes three different cases: Scalar pentagons (i.e. scalar five-point 
functions), scalar hexagons (scalar six-point functions) and scalar integrals with more than six 
propagators. 

Let us start with the pentagon. A five-point function in D = 4 — 2e dimensions can be expressed 
as a sum of four-point functions, where one propagator is removed, plus a five-point function in 
6 — 2e dimensions [102]. Since the (6 — 2e) -dimensional pentagon is finite and comes with an 
extra factor of £ in front, it does not contribute at 0(£"). In detail we have 

/5 = -leBlt^' - f b/J^ = - f + 0{e), (3.5) 

i=l i=l 

where I^^^^ denotes the (6 — 2£) -dimensional pentagon and denotes the four-point function, 
which is obtained from the pentagon by removing propagator /. The coefficients B and bj are 
obtained from the kinematical matrix Stj as follows: 

The six-point function can be expressed as a sum of five-point functions [103] without any correc- 
tion of 0(e) 

k = -tb4\ b, = l^{S-')^., (3.7) 
'■=1 j 

where the coefficients bi are again related to the kinematical matrix Sij. For the seven-point function 
and beyond we can again express the «-point function as a sum over {n — l)-point functions [106]: 



= -l^rilHl,. (3.8) 

;=1 



In contrast to eq. (3.7), the decomposition in eq. (3.8) is no longer unique. A possible set of 
coefficients r, can be obtained from the singular value decomposition of the Gram matrix 

4 

Gij = Y,UikWk{v'') (3.9) 



as follows [107] 



ri = ^, l</<«-l, r„ = -£o-, W5 = \Y^GjjVj5. (3.10) 
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3.1.2 Improvements of the Passarino-Veltman algorithm 

We now consider the reduction of tensor loop integrals (e.g. integrals, where the loop momen- 
tum appears in the numerator) to a set of scalar loop integrals (e.g. integrals, where the numerator 
is independent of the loop momentum). For one-loop integrals a systematic algorithm has been 
first worked out by Passarino and Veltman [108]. Consider the following three-point integral 

k2{k-piy{k-py-p2)^' ^ ^^^^ 

where p\ and p2 denote the external momenta. The reduction technique according to Passarino and 
Veltman consists in writing /^^ in the most general form in terms of form factors times external 
momenta and/or the metric tensor. In our example above we would write 

73^^^ = pi'prC21+p^P2'C22 + K,;^2'}C23+g^'C24, (3.12) 

where {Pi,P2} = p'iP2 + P2Pi- ^"^^ ^^^""^ solves for the form factors C2i,C22,C23 and C24 by first 
contracting both sides with the external momenta PiP\, P2P2, {Pi ,P2} and the metric tensor g^^^. 
On the left-hand side the resulting scalar products between the loop momentum k^ and the external 
momenta are rewritten in terms of the propagators, as for example 

2pi-k = k^-{k-pif + pl (3.13) 

The first two terms of the right-hand side above cancel propagators, whereas the last term does not 
involve the loop momentum anymore. The remaining step is to solve for the form-factors C2, by 



inverting the matrix which one obtains on the right-hand side of equation ( |3.12[ ). Due to this step 
Gram determinants usually appear in the denominator of the final expression. In the example above 
we would encounter the Gram determinant of the triangle 



pI Pi -Pi 
Pi- Pi p\ 



(3.14) 



One drawback of this algorithm is closely related to these determinants : In a phase space region 
where p\ becomes collinear to p2, the Gram determinant will tend to zero, and the form factors 
will take large values, with possible large cancellations among them. This makes it difficult to 
set up a stable numerical program for automated evaluation of tensor loop integrals. Quite some 
effort went therefore into improvements and alternatives, which avoid these instabilities. They are 
centered around the following ideas: 

• The use of a different reduction scheme in critical regions, based on an expansion around the 
small invariants [109-112]. 

• Choosing a different set of basic integrals. The new set can be non-minimal. Quite often the 
new basic integrals con^espond to scalar integrals in higher dimensions. [107, 113-118]. 

• Direct numerical integration [119-124]. 

• The use of spinor techniques, which avoid to a certain extent the occurrence of Gram deter- 
minants [125-130]. 
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Figure 7: Double, triple and quadruple cuts. 

3.1.3 Methods based on unitarity 

We know that the final answer for a one-loop amplitude in a massless gauge theory can be 
written as 

Ai^) = Ic,-/« +£c,/(''^') + £ c,,-,/f ^) +R. (3.15) 

; (,;■ i.j.k 

h, h and ^ are the scalar bubble, triangle and box integral functions. In a massive theory we 
would have in addition also scalar one-point functions. In a massless theory these functions are 
zero within dimensional regularisation. Note that there are no integral functions with more than 
four internal propagators. These higher-point functions can always be reduced to the set above, 
as we have seen in section [3.1.1[ R is called the rational term. The set of all occurring integral 
functions 

^ = (3.16) 

is rather easily obtained from pinching in all possible ways internal propagators in all occurring 
diagrams. We can assume that we know this set in advance. To compute the amplitude requires 



therefore the determination of the coefficients c,, c,j, cij^u and the rational term R. In eq. ( |3.15[ ) 
any £-dependence of the coefficients has been removed. The original £-dependent parts of the 
coefficients, which are multiplied with poles from the integral functions are collected in the rational 
term R. Within the unitarity-based methods one calculates the coefficients and the rational term 
without resorting to a Feynman diagram calculation. 

The original formulation of the method [131, 132] is based on the observation that the basic 
integral functions contain logarithms and dilogarithms, which develop imaginary parts in certain 
regions of phase space, for example 

Imlnf ^^^1 = -K[d{s)-e{t)\, 



^ -t-iO 

ImLiil l--r— — ?T ) = -lnfl--Vmln( ^^^^ ) . (3.17) 



{-t-iO) J \ tJ \-t-iO 

Knowing the imaginary parts, one can reconstruct uniquely the corresponding integral functions. 
In general there will be imaginary parts corresponding to different channels (e.g. to the different 
possibilities to cut a one-loop diagram into two parts). The imaginary part in one channel of a 
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one-loop amplitude can be obtained via unitarity from a phase space integral over two tree-level 
amplitudes. With the help of the Cutkosky rules we have 

ImAd) = Im / -f^^J^AM^. (3.18) 

A^') is the one-loop amplitude under consideration, A^*^^ and A^*^^ are tree-level amplitudes appear- 
ing on the left and right side of the cut in a given channel, as shown in the first picture of fig. ^. 



Lifting eq. ( 3.18 ) one obtains 



where "cut free pieces" denote contributions which do not develop an imaginary part in this partic- 
ular channel. By evaluating the cut, one determines the coefficients of the integral functions, which 
have an imaginary part in this channel. Iterating over all possible cuts, one finds all coefficients. 
One advantage of a cut-based calculation is that one starts with tree amplitudes on both sides of 
the cut, which are already sums of Feynman diagrams. Therefore cancellations and simplifica- 
tions, which usually occur between various diagrams, can already be performed before we start 
the calculation of the loop amplitude. This technique was used in the calculation of the one-loop 
amplitudes for e^e^ —>■ 4 partons [133, 134]. The rational part R can be obtained by calculating 
higher order terms in e within the cut-based method. At one-loop order an arbitrary scale /i^^ is in- 
troduced in order to keep the coupling dimensionless. In a massless theory the factor jx^^ is always 
accompanied by some kinematical invariant s^^ for dimensional reasons. If we write symbolically 

A'- = U f S V% ^c: (^y\J ^) ' , (3.20) 



the cut-free pieces co(5o/M^) ^ can be detected at order e: 



co(^) =co-ecoln( (3.21) 



With the advent of the new methods based on twistors, significant improvements were added to 
the unitarity-based technique [135-152]. Apart from the two-particle cut discussed above, one can 
also consider triple or quadruple cuts as shown in fig. ^ A particular nice result follows from 
quadruple cuts [135]: The coefficients of the box integral functions are given as a product of four 
tree amplitudes, summed over the two solutions of the on-shell conditions 

ll = ll = ll = ll= 0. (3.22) 

Unitarity-based methods contributed significantly to the calculation of the one-loop six-gluon am- 
plitude [131, 132, 136-140, 143, 145-147, 153-157] and the one-loop six-photon amphtude [152, 
158-161]. 
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3.2 Two-loop amplitudes and beyond 

For NNLO calculation two-loop amplitudes are required. The relevant two-loop integrals 
are far from trivial. In this section I review a few techniques which can used for the calculation 



of two-loop integrals. The Mellin-Bames transformation is discussed in section 3.2. 1|. Multiple 



poly logarithms together with their algebraic properties are introduced in section 3.2.2. Section 



3.2.3| is devoted to sector decomposition. Not treated in detail, but equally important are methods 



which reduce the work-load: Integration-by-parts identities [162] have a long-standing tradition. 
In addition there are the reduction algorithms of Tarasov [163, 164] and Laporta [165], which 
can be used to relate all tensor integrals to a set of basic scalar integrals. With the help of these 
techniques many two-loop amplitudes have been calculated: Bhabha scattering [166], pp —>■ 2 jets 
[167-173], e^e^ —>■ 3 jets [174-176] and Higgs production [177, 178]. In addition, the three-loop 
splitting functions, required for the evolution of parton distribution functions at NNLO, have been 
calculated [179, 180]. 

3.2.1 The Mellin-Barnes transformation 

The Mellin-Barnes transformation can be used to transform any Feynman parameter integral 
into a particular simple form, such that the integrals over the Feynman parameters can be per- 
formed: 

n r(vy) 

^ ^ (3.23) 






The Mellin-Barnes transformation reads 



r(vi + ... + v„)' 



(Ai+A2 + ...+A„) ' = ^—-—^ I doi... I dOn-i (3.24) 



1 

(271/) 

xr(-ai)...r(-a„_i)r(ai + ... + a„_i+c)Af...A,';r;A,7' 

Each contour is such that the poles of r(— a) are to the right and the poles of T{o + c) are to 
the left. This transformation can be used to convert any Feynman parameter integral to the form 



— /OO — /OO 

iCl 4 CT/i-l 4 -0-i-...-fT„_i- 



of eq. (3.23). Therefore we exchange the Feynman parameter integrals against multiple complex 
contour integrals. As this transformation converts sums into products it is the "inverse" of Feynman 
parametrisation. The contour integrals are then performed by closing the contour at infinity and 
summing up all residues which lie inside the contour. Here it is useful to know the residues of the 
Gamma function: 

res (r(a + a),a = -a-«) = ^ — res (r(-a + a),a = a + «) = — ^. (3.25) 

n! n'. 

Therefore we obtain (multiple) sum over residues. Techniques to manipulate these sums are dis- 
cussed in the next section. In particular simple cases the contour integrals can be performed in 
closed form with the help of two lemmas of Barnes. Barnes first lemma states that 



/OO 

doT{a + o)T{b + o)T{c-o)T{d -o) 



T{a + c)T{a + d)T{b + c)T{b + d) 



2ni J ^ ^ ^ ^ ^ ^ ^ ' r(a + Z7 + c + <i) ' 

(3.26) 



-/OO 
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if none of the poles of r(a + (J)r{b + <7) coincides with the ones from r(c — a)r{d — g). Barnes 
second lemma reads 



iTTi J 



^^ T{a + o)T{b + a)r(c + o)T{d - o)T{e - o) 



T{a + b + c-\-d + e-\-G) 

— loa 

r{a + d)T{b + d)T{c + d)T{a + e)T{b + e)T{c + e) 



(3.27) 



T{a-\-b-\-d + e)T{a-\-c + d + e)T{b-\-c + d-\-e) 

Although the MelUn-Bames transformation has been known for a long time, the method has seen 
a revival in applications in recent years [181-196]. 

3.2.2 Multiple polylogarithms 

The multiple polylogarithms are defined by 

il>!2>...>!t>0 '1 

The multiple polylogarithms are generaUsations of the classical polylogarithms Li„ (x) whose most 
prominent examples are 

Lii(x) = £ — = - ta(l -;c), U2{x) =Y.^, (3.29) 

;i=l '1 (1=1 h 

as well as Nielsen's generaUsed polylogarithms 

Sn,p{^) = Li„+i,i,...,i(ji:,l^^^), (3.30) 

p-i 

and the harmonic polylogarithms 

Hm\,...,mk^X) = Linij^,,,^/K|t(-'^; 1; — 1 1^ )- (3.31) 

k-l 

Multiple polylogarithms and the closely related harmonic sums have been studied extensively in 
the literature [197-214]. 

In addition, multiple polylogarithms have an integral representation. To discuss the integral 
representation it is convenient to introduce for 7^ the following functions 

G{zu...z,,y) = f-^ 7-^. (3.32) 

J ti-ZlJ t2-Z2 J h-Zk 


In this definition one variable is redundant due to the following scaling relation: 

G{zi,...,Zk;y) = G{xzi,...,xzk;xy) (3.33) 

To relate the multiple polylogarithms to the functions G it is convenient to introduce the following 
short-hand notation: 

Gmu...,mk{zi,-,zk;y) = G{0^^^^,zi,...,zk-i,0^-^,zk;y) (3.34) 

mi — I wj— 1 
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Figure 8: Sketch of the proof for the quasi-shuffle product of nested sums. The sum over the square is 
replaced by the sum over the three regions on the rh.s. 




Figure 9: Sketch of the proof for the shuffle product of two iterated integrals. The integral over the square 
is replaced by two integrals over the upper and lower triangle. 



Here, all zy for 7 = 1, are assumed to be non-zero. One then finds 



Limi mi,{xi,...,Xk) = (— 1)^G,„, mi, ( — , ;1 ] • (3.35) 

* Xi X1X2 X\...Xk 



Eq. ( [3.35[ ) together with ( [3.34| ) and ( [3.32| ) defines an integral representation for the multiple poly- 
logarithms. Multiple polylogarithms form two (independent) algebras. As a consequence, products 
of multiple polylogarithms can again be expressed as sums of single polylogarithms. To give an 
example, one has 

Lijni (-^1 )Lim2 ('"-2) — Li^j ,^2 (-"-l 1X2) -\- \A,fi2ffn {x2,X\ ) + Limi +HI2 ix\X2')i 

G{zuy)G{z2\y) = G{z\,Z2\y) + G{z2,z\;y). (3.36) 

The first relation is based on the representation in terms of nested sums. A sketch of the proof is 
shown in fig. The second relation is based on the integral representation. The corresponding 
sketch of the proof is shown in fig. ^. These algebraic properties allow automated manipulations of 
expressions related to multiple polylogarithms by computer algebra programs. Several programs 
are available [202,215-218]. 

Up to now we treated multiple polylogarithms from an algebraic point of view. Equally im- 
portant are the analytical properties, which are needed for an efficient numerical evaluation. As an 
example I first discuss the numerical evaluation of the dilogarithm [219]: 



X 



r \ni\-t) ^ x" 
Li2(x) = - / dt \ > = £ - (3.37) 
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Figure 10: Sector decomposition splits the integration over the large triangle xi +X2 < 1 into the sectors 
xi > X2 and X2 > xi . 



The power series expansion can be evaluated numerically, provided |x| < 1. Using the functional 
equations 

Li2(x) = -Li2 0)-^-l(ln(-x))^ 

Li2(x) = -Li2(l-x) + — -ln(x)ln(l-x). (3.38) 

6 

any argument of the dilogarithm can be mapped into the region |x| < 1 and — 1 < Re(x) < 1 /2. The 
numerical computation can be accelerated by using an expansion in [— ln(l — x)] and the Bernoulli 
numbers B, : 

OO T> 

,=0 '-J- 

The generalisation for the numerical evaluation of multiple polylogarithms has been worked out 
in [211]. 



3.2.3 Sector decomposition 

In this section I will discuss an algorithm, which allows to compute numerically the coeffi- 
cients of the Laurent expansion for a multi-loop integral for a given kinematical configuration of 
external momenta. The major challenge such an algorithm has to face is the disentanglement of 
overlapping singularities. An example for an overlapping singularity is given by 



(3.40) 



The term 1 / (xi +X2) is an overlapping singularity. Sector decomposition [220-222] is a convenient 
tool to disentangle overlapping singularities. For the example in eq. ( p. 40 ) one splits the integration 




region into two sectors xi > X2 and xi < X2. This is shown in fig. 10 In the first sector one 
rescales X2 as x'2 = X2/X1, while in the second sector one rescales x[ = x\jx2. By applying the 
sector decomposition iteratively, one finally arrives at a form where all singularities are factorised 
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Figure 11: Cancellation of infrared divergences between virtual corrections and real corrections. 



explicitly in terms of factors of Feynman parameters like Xj ^. Subtractions of the form 

1 1 

J dxjxj'-'f{xj) = -1 /(O) + J dxjxj'-' [/(xy) -/(O)] (3.41) 


for each j, where limxj^of{xj) is finite by construction, allow to extract all poles and lead to 
integrals which are finite and can be integrated numerically. 



4. Cancellation of infrared divergences 

Infrared divergences occur at next-to-leading order and beyond. At NLO real and virtual cor- 
rections contribute. The virtual corrections contain the loop integrals and can have, in addition 
to ultraviolet divergences, infrared divergences. If loop amplitudes are calculated in dimensional 
regularisation, the IR divergences manifest themselves as explicit poles in the dimensional regu- 
larisation parameter e = 2 — D/2. These poles cancel with similar poles arising from amplitudes 
with additional partons but less internal loops, when integrated over phase space regions where 



two (or more) partons become "close" to each other. This is illustrated in fig. |11| In general, the 
Kinoshita-Lee-Nauenberg theorem [223, 224] guarantees that any infrared-safe observable, when 
summed over all states degenerate according to some resolution criteria, will be finite. However, the 
cancellation occurs only after the integration over the unresolved phase space has been performed 
and prevents thus a naive Monte Carlo approach for a fully exclusive calculation. It is therefore 
necessary to cancel first analytically all infrared divergences and to use Monte Carlo methods only 
after this step has been performed. 

4.1 Infrared divergences at NLO 

At NLO, general methods to circumvent this problem are known. This is possible due to the 
universality of the singular behaviour of the amplitudes in soft and collinear limits. Examples are 
the phase-space slicing method [225-227] and the subtraction method [228-232]. It is worth to 
examine a simple NLO example in detail to understand the basic concepts. We consider the NLO 
corrections to 7* ^ 2 jets. The real corrections are given by the matrix element for 7* qgq and 
read, up to colour and coupling factors 



= 8(1 



2 2 2,, ,X2 -Xi 

+ (l-£)— + (l-£)— -2£ 

X1X2 X\ X2 X\ X2 



(4.1) 
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where xi = >si2Ai23 andx2 = J23Ai23- This term is integrated over the three particle phase space. 
Singularities occur at the boundaries of the integration region at xi =0 and X2 = 0. Within the 
subtraction method one subtracts a suitable approximation term da^ from the real corrections 
da^. This approximation term must have the same singularity structure as the real corrections. If 
in addition the approximation term is simple enough, such that it can be integrated analytically over 
a one-parton subspace, then the result can be added back to the virtual corrections da^. 



^NLO 



- J dG'' + JdG''= I [da^-dO^)+ I [jC7^ + JdoA 
n+l n n+l n \ 1 / 



(4.2) 



Since by definition da^ has the same singular behaviour as da^, da^ acts as a local counter-term 
and the combination {da^ — dc^) is integrable and can be evaluated numerically. Secondly, the 
analytic integration of do^ over the one-parton subspace will yield the explicit poles in e needed to 
cancel the corresponding poles in do"^ . For the example discussed above the approximation term 
can be taken as a sum of two dipole subtraction terms: 



+ \s^M,p'i)\' ^ 



Xi(xi 4-X2 

2 



Xl X\ 



(4.3) 



The momenta p\, p'^, p'[ and are linear combinations of the original momenta p\, pi and p^. 
The first term is an approximation for xi 0, whereas the second term is an approximation for 
X2 0. Note that the soft singularity is shared between the two dipole terms and that in general the 
Bom amplitudes are evaluated with different momenta. Antenna subtraction [233,234] allows 
to reduce the number of subtraction terms needed and interpolates smoothly between the xi ^ 
and X2 regions. Within antenna subtraction one could take for our example as approximation 

"222 

Xl 



k2(M",K')r 



1 



^123 



^X\ 



XlX2 



X2 



Xl ' X2 



(4.4) 



Again, the Born amplitude is evaluated with momenta p'y and p'^, which are linear combina- 
tions of the original momenta pi, pi and ps. As a rule of thumb, one antenna subtraction term 
equals two dipol subtraction terms. 

The dipole subtraction scheme has been worked out for general NLO calculations. The matrix 
element corresponding to the approximation term da^ is given as a sum over dipoles: 



pairs ij kj^ij 



(4.5) 



Each dipole contribution has the following form: 

^.._^(0)*/ ~. ~ N (-Tfc- T,-y) Vij,k ^(0) 

■^ij,k — \Pl->---^P{i])^---^Pk,---) „2 



Q+'2{pi,.:,P(ij),-,Pk,--) ■ (4.6) 



Here T, denotes the colour charge operator for parton i and Vij^k is a matrix in the spin space of the 
emitter parton (ij). In general, the operators T,- lead to colour correlations, while the Vij/s lead to 
spin correlations. The generation of all subtraction terms can be automated [81]. 
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4.2 Infrared divergences at NNLO 

The following terms contribute at NNLO: 

doP = * s^P + s^P * ^P + ^i^) * ^i^)) (4.7) 

where ^/n '^ denotes an amplitude with n external partons and / loops, d^n is the phase space 
measure for n partons. We would like to construct a numerical program for an arbitrary infrared 
safe observable O. Several options for the cancellation of infrared divergences have been discussed 
[235-253]. Among those, the subtraction method - well-known from NLO computations and 
sector decomposition are the most promising candidates. Let us look again at the subtraction 
method. To render the individual contributions finite, one adds and subtracts suitable pieces: 

+ J On+i dal^^^ + On+iodaf;'^ -Onodai'^'^ 
+ jOn d(yP + 0„odaj!^'^^ + 0„odaj,^'^\ 

Here daj^^\^ is a subtraction term for single unresolved configurations of Born amplitudes. This 
term is already known from NLO calculations. The term dajj^'^^ is a subtraction term for double un- 
resolved configurations. Finally, daji^'^^ is a subtraction term for single unresolved configurations 
involving one-loop amplitudes. 

To construct these terms the universal factorisation properties of QCD ampUtudes in unre- 
solved limits are essential. QCD amplitudes factorise if they are decomposed into primitive ampli- 
tudes. Primitive amplitudes are defined by a fixed cyclic ordering of the QCD partons, a definite 
routing of the external fermion lines through the diagram and the particle content circulating in the 
loop. One-loop amplitudes factorise in single unresolved limits as [131,236,254-258] 

A« = Sing(O'i) .^(1)^ +sing(i'i) -A^^l,. (4.8) 

Tree amplitudes factorise in the double unresolved limits as [235, 259-264] 

Ai«) = Sing(«.2) .Af_V (4.9) 

In addition, the pole structure of two-loop amplitudes is known [265-267]. The subtraction terms 
have to interpolate between the various singular limits. Spin-averaged antenna subtraction terms 
have been worked out in [60, 245] and first results at NNLO based on the subtraction method are 
available for e~^e~ 2 jets [60,61] and the thrust distribution [62]. 

5. Summary 

Precision calculations for multi-parton processes will play an important role for the physics 
program at the LHC. I discussed various approaches how to overcome the bottle-necks we face 
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technically: the length of expressions, the calculation of loop integrals and the occurrence of infra- 
red divergences. Computer algebra, a divide-and-conquer approach and new developments based 
on on-shell recursion relations help us to keep the length of expressions manageable. Sophisti- 
cated techniques related to Mellin-Barnes transformations, shuffle algebras or sector decomposi- 
tion allow us to compute multi-loop integrals. Urgently needed are also automated computations 
of one-loop amplitudes. I discussed improvements of the Passarino-Veltman algorithm and new 
ideas base on generalised unitarity. As far as the cancellation of infrared divergences is concerned, 
there are systematic methods at NLO. Several methods, which address the issue at NNLO have 
been proposed and first numerical results for observables based on the use of these methods have 
emerged. 
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